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ABSTRACT 


The Burrell Orbital Flux and Energy Spectra Code calculates the position 
of an earth satellite in geocentric and B-L coordinates. Utilizing the B-L co- 
ordinates, the code calculates the radiation flux and energy spectra for protons 
or electrons at any point in a circular or elliptical orbit. 
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TECHNICAL MEMORANDUM X-53406 


ORBITAL CALCULATIONS AND TRAPPED RADIATION MAPPING ‘ 

SUMMARY 


Given the orbital parameters of an orbiting satellite, the BOFES code 
calculates the position of the satellite in geocentric and B-L coordinates. 
Satellite positional data have been generated by the code duplicating the ephem- 
eris data for the Pegasus satellites for a period of two weeks. 

The code allows the complete mapping of the geomagnetic field traversed 
by the satellite on any feasible elliptical orbit out as far as the lunar orbit. 

The code calculates the radiation intensity and energy spectral of trapped 
radiation at each orbital step and obtains a time flux weighted average at the 
end of each orbit. Later, a method will be added to the proton and also to the 
electron code to calculate the dose rates behind several shield thicknesses at 
each step in the orbit. 


INTRODUCTION 


The Burrell Orbital Flux and Energy ^ectra (BOFES) code gives an 
accurate calculation of the orbital position of earth satellites and simultaneously 
computes the geomagnetic B-L coordinate of the satellite. The B-L coordinates 
are calculated using the 48-term representation of the geomagnetic field by 
Jensen and Cain or a simple dipole representation. The BOFES code utilizes 
the B-L coordinates to determine the radiation fluxes and energy spectra at any 
point of the orbit. The code allows the complete mapping of the trapped radia- 
tion field traversed by the spacecraft on any feasible elliptical or circular orbit 
around the earth and out as far as the lunar radius. The number of orbits and 
also the number of steps per orbit are chosen at the convenience of the users. 


COORDINATE TRANSFORMATION 


To determine the position of a satellite, one must be given a set of six 
orbital elements that specify the orbit's orientation in space, its size and shape. 



and the position occupied by the object at a specified time or the time at which 
it is at a specified point. For geocentric Orbits, such a set consists of three 
orientation elements and three dimensional elements. 

The orbital elements are 

(1) Q, the longitude of the node, measured in the plane of the equator 
from the direction of the vernal equinox eastward to the direction of the ascend- 
ing node, or intersection of the orbit with the equator. 

(2) cu, the argument of perigee, the angle measured in the plane of the 
orbit from the direction of the ascending node to the direction of perigee. 

(3) 0, the Inclination, the angle between the equatorial plane and the 
plane of the orbit. 

The three remaining "dimensional" elements are 

(1) a, the semimajor axis 

(2) e, the eccentricity, the ratio of the distance from the center of the 
orbit to the focus (center of earth) to the semimajor axis. 

(3) T , the time of perigee passage. 

P 

The orientation elements, n, <u, <p, are typical Euler angles. Using 
these angles a coordinate transformation can be made from tiie orbital system 
to the earth coordinate system, to obtain the latitude and longitude. Using the 
Euler angles, the elements a^j of the coordinate transformation matrix are 
givqn by 


a^i = 
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sin u) 

sin £2 cos 0 

^12 “ 
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■ cos ct! sin £2 cos 0 
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a 23 = -sin (p cos n 


a 3 j = sin o) sin <p 

^32 ~ ^ 0 

a 33 ~ CO s (p 


(1) 


The coordinate transformation to the earth system (x, y, z) is obtained 
by the following matrix multiplication: 


y 

j 


/an 

ai2 

ai3 

a2i 

a22 

a23 

\a3i 

as2 

as3 

is the coordinate 



( 2 ) 


drawn from the pericenter of the ellipse to the perigee and z* is the axis 
normal to the plane of the ellipse. The z* element in the matrix of equation (2) 
is defined as zero since the point on the orbit is in the plane of the ellipse. 

From equation (2) the geocentric coordinate system is given by 

x= aji X* + aj 2 y* 

y = a 2 i X* + a 22 y* (3) 

z - a 3 i X* + a32 y* 

Using the geometry of Figure 1 and spherical coordinates, the point P (x* , y=!= ) 
is placed in the (x, y, z) system fcy the following: 

X = r sin X' cos ip' 

y = r sin X' sin ip' (4) 

z = r cos X' 

Substituting the above into equation (3) we obtain: 

. z ^31 + ^32 y* 

cosX' = -= . (5) 

r r ' ' 
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VERNAL EQUINOX 


FIGURE 1. COORDINATE TRANSFORMATION 
Since X is the latitude of point P (x* , y* ) measured from the equatorial plane, 
then 

\' = 90 - A 


and 

cos A' = cos (90 -A) = sin A. 
Thus equation (5) becomes 

^31 ^ ^32 y* 

sin A = . 

r 


( 6 ) 
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From Figure 1 and equation (3) 


tan ^ 

^ X 


\ ^ ^2 y* 

an X* + an y* ’ 


(7) 


where tj)' is a displaced longitude (in the equatorial plane) of the point P (x* , y*) 
measured eastward from the direction of the vernal equinox to the line of inter- 
section of the projection of the Vector r onto the equatorial plane. Thus equations 
(6) and (7) give the latitude and a displaced longitude in the equatorial plane 
measured from the vernal equinox. To find the true longitude, one must know 
the angle "o!" between the Greenwich Meridian and the vernal equinox at a spe- 
cific time. 

Thus the true longitude ip is given by 

Ip = ip' - a . (7a) 

The method of determining a is discussed in the next section on methods. 


METHOD 


The BOFES code utilizes a constant angle step during an orbit. At each 
step the position relative to the earth coordinate system is computed. Simul- 
taneously, the B-L coordinate of each point is also calculated. The B-L coordi- 
nate is then used to calculate the radiation flux and energy spectra from Input 
data. The details of the methods used follow. 

Given the semimajor axis a the period of the orbit is determined by 

T = 27t \/^ ; p = GM, (8) 

where G is the gravitational constant and M is the mass of the earth. 

The radial distance of the satellite (Fig. 2) and the flight time (At) for 
each constant step in angle (A©) are given by, ©k+j = ©!{■*■ where 0 q = 0° 
at perigee 
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( 9 ) 


-g - 


fc<-l i + e cos 0 


k+1 


a - r 


y = 


k+1 


a e 


[cos-i(y) - cN/rr" 


( 10 ) 


. gg 


At = 


Vl 



where is the radial distance measured from pericenter to the position of 
the spacecraft, Gj^+j is the angle measured from perigee to the spacecraft 
position in the orbit plane, e is the eccentricity of the orbit, y is equivalent to 
the cosine of the eccentric anomaly and tj,+ j is the fli^t time of the spacecraft! 
corresponding to the angle Using equation (6) , the radius of the earth for 

the latitude, A, is computed by 

0.996633 Re 

R— ^ 1 

n/i - 0. 0067227 cos^ A 



where Re is the equatorial radius [ 1] . The altitude of the spacecraft above the 
earth’s surface is given by 


H = r - R, 


(13) 



where r is obtained from equation ( 9) ; units of the altitude are the same as Re 
and the semimajor axis a. 

The geodetic latitude is computed [ 1] , using equation ( 6) , as 


X , = tan“^ 
d 


" tan X 

0. 9932773 ‘ 

C -J 


(14) 


The longitude il>, measured eastward from the Greenwich Meridian 
(Fig. 3) , is given by 


^ = il>' - a 


(0 ip < 2ir). 


(15) 


NORTH 



The angle o; (in radians) is given by 

a = ib + 6. 300388 t, (16) 

o 

where the constant is the rate of change of the mean sidereal time in radians/ day; 
t is the time of flight (in days) measured from beginning of the calculations; and 
is angular distance (in radians) between the vernal equinox and the Greenwich 
Meridian at the beginning time of the calculation. The ip^ is given by 
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(17) 


tb = . 98564995 [ MJD - 35 838. 0] + 99. 39235, 

where MJD is the mean Julian date of the beginning of the calculation; the con- 
stant 99. 39235 is >Po at MJD 35 838. 0; and 0. 98564995 is the rate of change of the 
mean sidereal time in degrees/day. (The times of epoch are given in Julian 
days, and for convenience 2 400 000. 5 has been subtracted to give a "modified" 
Julian day. ) From equations ( 13) , ( 14) , and (15) , the position of the space- 
craft is determined at each point in time. 

In the BOFES codes , all computations are begun at the perigee point in 
the orbit. When given the orbital parameters for an ep>och Tq, which may or 
may not be the time of perigee passage, the time of perigee passage must be 
calculated and the argument of perigee o) and the longitude of the ascending 
node n must be corrected. The time of the previous perigee passage correspond- 
ing to the epoch Tq is given by 


TMq 


where To is the epoch time; Mq is the mean anomaly corresponding to Tq; and T 
is the period of the orbit. When Tp is computed, the orbital parameters co and 
S2 are computed for the time of perigee passage Tp. These computations are 

ii^ut data designated as 
date (days) , midnight 
Greenwhich time, at which computations are to be made and TIME is the frac- 
tion of a day (day units) after midnight when computations are performed. 

Using the position elements of the satellite, the B-L coordinates can be 
determined simultaneously ly utilizing the 48-term representation of the geo- 
magnetic field by Jensen and Cain or a simple dipole representation [2] may be 
obtained by using^ equations (18) and (19). Hence, 


performed by hand. The time of perigee passage T is 
the sum of DATE and TIME, where DATE is the JuUan 


L = 


r 

cos^ A 


and 






1/2 


(18) 


(19) 
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where A is the geocentric latitude given in equation (6) ; r is the radial distance 
(in earth radii) obtained by equation ( 9) ; and M is the dipole moment (. 311653 
gauss Re®). 

Since the radiation (proton and electron) fluxes are functions of the B-L 
coordinates, the BOFES code utilizes the B-L coordinates to determine the 
radiation intensity at any point in an orbit and as far out as the lunar radius. 

The code is also used to determine the artificial belt radiation level at any 
point on the orbit. 


PROTON CODE 


The proton BOFES code computes the proton flux 4>[ B(t) , L(t)] at each 
orbital step. Also a time average flux is determined continuously and printed 
at the end of each orbit. The time average flux is given by 
N 

2 $ [B(t). L(t)] At 

- i = 1 

( 20 ) 


<p - 


N 


Z 


i = 1 


N 

where ^ At is the accumulative time at the end of the step in the orbit, 
i = 1 

The B-L coordinate is also used to determine the corresponding exponential 
energy parameter Eq which is employed to find the differential and integral pro- 
ton energy spectrum at each point in the orbit calculation. 

After the proton flux, has been determined, the time-weighted proton 
differential and integral energy spectra are found for certain specified input pro- 
ton energies, E(n) . The code computes these energy spectra for each point in 
the orbit and keeps running totals for each proton energy, E(n) . The time- 
weighted average of the integral and differential energy spectra for the 
energy group are 


I (> E ) = 
n' 


N E - E(n) 


I 


( 21 ) 


e 



N (B,L) Ej^-E(n) 

$ (E ) = 

' n ^ V--1 


The above are printed out after each orbit for each proton energy E(n) . E]y[ is 
the minimum energy of the proton data, Eq is the exponential energy parameter 
corresponding to each point in the orbit, Atj is the time step along the orbit and 
the summation is the accumulation of the total flight time of the spacecraft. 


ELECTRON CODE 

The BO FES code for electrons computes the instantaneous electron flux 
at any point in the orbit and also carries a time averse of the electron flux 
given by equation (20) . 

The electron’ s energy spectra above 0. 5 MeV are given in eight energy 
bands and are normalized so that 

2 f (E) = 1. (22) 

j=l ^ 


The electron energy spectra are a function only of the L coordinate. Therefore, 
for any point in the orbit, the energy bands corresponding to the L coordinate are 
obtained, and an accumulative time-flux-weighted average is printed out at the 
end of each orbit and is given by 

^ $j(B,L) f.[E(t), L{t)l At. 

T.(E)=^^ . (23) 

J L)At. 


where 


j= 1. 2, ... 8. 

The energy groups represented above are the same as those in Reference 3, 
page 48 (i. e. , 0, 5-1, 1-2, 2-3 MeV, . . . ) . 
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SPECIAL FEATURES OF CODE 


The option MT determines the position of the orbit angle step relative to 
the time step. For MT =0, 1 the angle step along an orbit is halfway between 
the time step, or the angle step does not coincide with the time step. For 
MT = -1, the angle and time step coincide, i. e. , t^+i is the flight time for the 

orbit angle 0k+-l* 

The option JOB distingviishes the tsTJe of calculations to be performed and 
the output format to be used. When JOB =1, 2, or 3 the output is at each 
orbital step; for JOB = 4, the output is at the end of each orbit only. When 
JOB = 1, the code computes and outputs the position (degrees) of the spacecraft 
in orbit (6), the altitude (H), latitude (2d), longitude (ip) , the time step (At) 
and the total flight time (TT) , respectively. When JOB = 2, the output is the 
same as for JOB = 1 with the B and L coordinates included. For JOB = 3, the 
output is the same as for JOB =2, but includes the radiation flux (particles per 
^ a running sum of the product of flux and the time step (no. of 
particles crossing unit area during the time step) , also the reciprocal of the 
exponential energy parameter Eq (proton code). The reciprocal is printed as 
zero if the energy spectra option is 1 (k= 1). For JOB = 4, the output includes 
the total flight time from beginning of calculations (TT) , the product of the flux 
and the time step (SFLUX) and the time averaged particle flux (SFLUX/TT) 
given by equation (20) . 

An option A L is also Included which determines the maximum value of the 
L coordinate to be calculated using the 48-term representation of the geomagnetic 
field. All L values (in earth radii units) equal to or greater than the value of 
AL will be computed using a dipole representation. For altitudes below 800 km 
and longitudes between 90 and 240 degrees east, a 50 percent reduction [ 4] in run 
time can be realized if one bypasses the calculation.of the geomagnetic field using 
the 48-term expansion by Jensen and Cain. If the position of the spacecr^t falls 
within the above region, a test is provided to lypass the multipole expansion of 
the geomagnetic field and computes the dipole field and sets the radiation fl^ 
equal to 1 particle per cm^ -sec. An option to bypass this geogiraphic test is 
provided because, if one is interested in the geomagnetic field distribution, the 
B-L coordinates must be computed at each position traversed by the spacecraft. 

The energy spectra option k provides for the calculation of energy spectra 
(protons and electrons) at each orbital position. When k= 2, the electron 
energy bands corresponding to the L coordinate are determined and printed after 
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each orbital step and a time-flux-weighted average is printed at the end of each 
orbit. For protons (k = 2) , the time-flux-weighted average of the differential 
and integral energy spectra (Eq. 21) corresponding to the proton energies E(n) 
are computed and printed after each orbit. When k = 1, the energy spectra for 
protons and the eight bands (electrons) are not given as output. 


RESULTS AND CONCLUSIONS 


Figures 4 through 7 show the different types of computations possib^ 
using the BOFES (proton and electron) codes. 

Figure 4 is a plot of the proton flux at an altitude of 1000 km computed 
along a circular polar orbit. Curve B shows the total proton flux per orbit 
while curve A is a running average of the total flux per orbit, which tends to 
stabilize after several orbits. Figure 5 shows the proton flux on a circular 
orbit, inclined 32° . 5, at an altitude of 286 kilometers. Again, curve B is the 
total flux per orbit, which shows a greater fluctuation than Figure 4. Curve A 
is the running average of the proton flux. Similar electron data can be obtained 
using the electron flux data AEl. 

Figure 6 shows the electron flux plotted as a function of altitude for four 
different inclinations using the electron data (AEl). 

Figure 7 shows the proton flux as a function of altitude uging the proton 
data (Apl) . 

Given the orbital parameters of an orbiting satellite, the BOFES code can 
compute the satellite positional data duplicating the satellite ephemerous data 
from more elaborate codes with only small errors in the position of the satellite 
for a period of two weeks. 

Subsequently, a method will be added to the proton, and possibly to the 
electron, code whereby the radiation dose intensity behind several shield thick- 
nesses will be computed at each orbital step and a time and flux weighted average 
will also be computed for each orbit. 

The code allows the complete mapping of the tr^ped radiation field 
traversed by the satellite on any feasible elliptical or circular orbit around the 
earth and out as far as the lunar radius. 
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ORBIT NUMBER 


FIGURE 4. PLOT OF THE PROTON 
FLUX (protons per cm^-hr) 
CALCULATED ON A CIRCULAR 
POLAR ORBIT AT AN ALTITUDE 
OF 1000 km 



FIGURE 5. PROTON FLUX (protons per 
cm^-hr) CALCULATED ON A CIRCULAR 
ORBIT AT AN ALTITUDE OF 286 km 
INCLINED 32?5 WITH THE EQUATOR. 



ALTITUDE, nmi 

FIGURE 6. ORBITAL INTEGRATIONS 
WITH ELECTRON DATA AEl, E > 0. 5 
MeV FOR CIRCULAR ORBITS (REF. 

3, p. 97) 



ALTITUDE (nmi) 

FIGURE 7. ORBITAL INTEGRATIONS 
USING PROTON DATA API, E > 34 
MeV (REF. 3, p. 159) 
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APPENDIX 


Key to Proton Flux - Energy Spectra Data 

The following is a description of the data computed by the BOFES proton 
code (Fig. A-1) . The sample data is for printout option equal to 3 (JOB = 3) , 
which is inclusive of JOB =1,2, and 4, The energy spectra are not included in 
the output of JOB = 1 and 2. 

1. Orbital parameters described previously in the section Data Input, 
except LA TO and LONGO which are the latitude and longitude of the Initial 
staJrtlng position, respectively. 

2. Minimum energy (in MeV units) of the proton data. 

3. THETA - angle (degrees) measured from perigee of orbit to position 
of spacecraft. 

ALT - Altitude of spacecraft (same units as semimajor axis). 

LAT - Latitude (degrees) measured from the equator. 

LONG - Longitude (degrees) measured eastward from Greenwich 
Meridian. 

B - B coordinate (gauss units). 

L - L coordinate (Earth radii units). 

TT - Total orbital flight time of spacecraft measured from the 

beginning of computations (same units as time factor, TF). 

FLUX - Proton Flux (particles cm~^ sec~*) at spacecraft position. 

SFLUX - Proton current (particles cm“^) , essentially the product 
of the proton flux and time between each orbital step (DT) 
with a time factor correction. 

1/EO - Reciprocal of energy parameter (MeV~^) 
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4. Total orbital flight time (TT), proton current (SFLUX) and the 
quotient of SFLUX and TT (equivalent to proton flux with time unit equal to 
time factor TF) , respectively. This printout is taken from the last printout 
for each orbit and is the output from JOB = 4. 

5. Time and flux -weighted average of the proton energy spectra data 
including proton energy (E) in MeV units, the differential proton flux at energy 
E and the integral proton flux for proton energies equal to and greater than 
energy E. 

Key to Electron Flux-Energy Spectra Data 

The following is a description of symbols used in the output of the 
BOFES Electron code (Fig. A -2). The sample data is for output option 3 
(JOB = 3) which includes the output of JOB =1, 2, or 4. The elctron energy 
spectra can be included only in the output from options 3 and 4. 

1. Same as 1 for proton data. 

2. Same as 3 for proton data. 

3. Eight differential energy spectra groups (. 5- IMeV, 1-2, 2-3,... 
6-7, 7-«>) , respectively, satisfying the condition of equation (22) . The dif- 
ferential energy spectra are a function of the L coordinate only, therefore, the 
value of the energy spectra in each of the eight groups is determined by the L 
value printed on the line above the energy spectra. 

4. Same as 4 for proton data, 

5. Time and flux-weighted average of the differential energy spectra 
data for each point along the orbit. 
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ORBITAL PARAMETERS FOR BOFES PROTON CODE 
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35.00 .524 + 03 .581 + 05 

50.00 .395 + 03 .445 + 05 

65.00 .303 + 03 .341 + 05 

80.00 .232 + 03 ,261 + 05 

100.00 .163 + 03 ,183 + 05 



ORBITAL PARAMETERS FOR BOFES ELECTRON CODE 
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Data Input - Electron Code 


The following table is a listing of the data input for the BOFES electron 

table. 


Variable 

Format 

Description 

(A) MAX 

15 

Number of L values for 
electron flux data. 

(B) KMAX (1), BL (1) 

I5,F7.0 

KMAX (1) is the number of 
B-flux pairs for L value, 
BL (1). 

(C) B(J,I), FL(J,I) 

F6. 0, E8. 0, 
3(F7, 0, E8. 0) 

B-flux pairs for L value , 
BL(1) card (B) is repeated 
until KMAX (max) BL(max) 
is read. 

(D) LMAX 

15 

Number of L values for 
energy spectra data. 

(E) BEL(l) 

F5. 2 

L coordinate value. 

(F) E (1.1) 
• 

8E10. 3 

Energy spectra, 8 energy 
bands . 5-1 MeV, 1-2 MeV, 
... 6-7 MeV, 7-«> MeV. 
card (E) is repeated. Until 
Bel (LMAX) is read. 

E (8,1) 



(G) NCASE 

15 

Number of different orbital 
cases to be rim. 

(H) A MEG, A MEGA, PHI 

3E14. 8 

A MEG is argument of perigee 
(degrees) A MEGA, longitude 
of ascending node (deg) , 
measured eastward from 
vernal equinox, PUl, incli- 
nation of orbit. ( degrees) . 
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Variable Format Description 

(I) DMEG, DMEGA, DATE, 

TIME 4E14. 8 DMEG, rate of change of 

AMEG (degrees/hr) 

DMEGA, rate of change of 
A MEGA (degrees/hr) . 

DATE, Julian date (days), 
midnight Greenwich time of 
day at which computations 
are to be made, TIME, 
fraction of day (day) after 
midnight when computations 
are performed. 

(J) A, ECC, TF, REE 4E14. 8 A is the semimajor axis of 

orbit ECC is the eccentricity 
TF is a time factor, gives 
time unit used in computa- 
tions. Ratio of (time unit 
used/ hr) 

TF = 60, time units in 
minutes 

TF = 1, time units in hours 
TF = 1/24, time unit in days, 
etc. 

REE is equatorial radius of 
earth (A and REE must have 
same units). 

(K) NC NS, ERR 2I5, F7.0 NC, No. of orbits per case 

(having same orbital param- 
eters). NS, no. of angle 
steps per orbit. ERR, ac- 
curacy of L computation eg. 
ERR = . 1 accuracy = 1% 

ERR = 1. accuracy = 

10% etc. 
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Variable Format Description 


(L) MT, JOB, AL. K 215, F7.0, IS MT, position of orbital angle 

relative to time step. JOB, 
specifies type of orbital cal- 
culations to be made. AL, 
upper limit on L computed 
by 48-term representation of 
geomagnetic field. K, 
Energy spectra option, 
previously described. 
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Data Input - Proton Code 


A listing of the input data for the BOFES proton code is given in the 
following table. 


Variable 

Format 

Description 

MAX 

15 

Same as (A) in electron code. 

KMAX (1) , BL(1) 

15, F7. 0 

Same as (B) in electron code. 

B(J,I), FL(J,I) 

F6. 0, E8. 0 
(F7. 0, E8. 0) 

Same as ( C) in electron code . 

EMIN 

FI 0.0 

Proton energy (MeV) of 
exponential parameter 
(EO) data. 

LMAX 

15 

Same as (D) in electron code. 

KEMAX (1), BEL(l) 

15, F7.0 

KEMAX, no. of B-EO pairs 
for L coordinate, BeL (1), 

BB(1,1), EO(l.l) 

F6.0, E8. 0 

B-exponential energy 
parameter (EO) pairs for L 
coordinate, BeL (1). 

BB(KEMAX,1), EO(KEMAX, 1) 


Repeated until 

KEMAX(LMAX) BeL(LMAX) 
is read. 

NCASE 

15 

Same as (G) in electron code. 

A MEG, A MEGA, PHI 

3E14. 8 

Same as (H) in electron code. 

DMEG, DMEGA, DATE, 
TIME 

4E14. 8 

Same as (I) in electron code. 

A, ECC, TF, PEE 

4E14.8 

Same as ( J) in electron code. 

NC, NS, ERR 

215, F7. 0 

Same as (K) in electron code. 

MT, JOB, AL, K 

215, F7. 0, 15 

Same as (L) in electron code. 
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